Effect of dipolar interactions on the magnetization of a cubic array of nanomagnets 
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We investigated the effect of intermolecular dipolar interactions on an ensemble of f 00 3D-systems 
of 5 x 5 x 4 nanomagnets, each with spin S = 5, arranged in a cubic lattice. We employed the Landau- 
Lifshitz-Gilbert equation to solve for the magnetization curves for several values of the damping con- 
stant, the induction sweep rate, the lattice constant, the temperature, and the magnetic anisotropy. 
We find that the smaller the damping constant, the stronger the maximum induction required to 
produce hysteresis. The shape of the hysteresis loops also depends on the damping constant. We 
find further that the system magnetizes and demagnetizes at decreasing magnetic field strengths 
with decreasing sweep rates, resulting in smaller hysteresis loops. Variations of the lattice constant 
within realistic values (f .5 nm - 2.5 nm) show that the dipolar interaction plays an important role 
in the magnetic hysteresis by controlling the relaxation process. The temperature dependencies of 
the damping constant and of the magnetization are presented and discussed with regard to recent 
experimental data on nanomagnets. Magnetic anisotropy enhances the size of the hysteresis loops 
for external fields parallel to the anisotropy axis, but decreases it for perpendicular external fields. 
Finally, we reproduce and test a previously reported magnetization curve for a 2D-system [M. Kayali 
and W. Saslow, Phys. Rev. B 70, f 74404 (2004)]. We show that its hysteretic behavior is only 
weakly dependent on the shape anisotropy field and the sweep rate, but depends sensitively upon 
the dipolar interactions. Although in 3D systems, dipole-dipole interactions generally diminish the 
hysteresis, in two-dimensional systems, they strongly enhance it. For both square two-dimensional 
and rectangular three-dimensional lattices with B\\(x + y), dipole-dipole interactions can cause large 
jumps in the magnetization. 

PACS numbers: 75.40.Mg, 75.60.Ej, 75.75.+a, 75.50.Xx 



INTRODUCTION 



The need of smaller memory storage devices, |1|, _ 
S&EHHUHmHElQthe interest in devel- 
oping quantum computing, |l5j and the hope for under- 
standing the relationship between the macroscopic and 
microscopic magnetic behaviors has led intense research 
into th e prop erties of nanoscale magnets. B El, E L T 
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Many issues still remain 
unclear and serious problems must be overcome in order 
for them to be technologically useful. Prominent among 
these is the loss of memory during magnetic relaxation. 

Ferromagnetic nanodots are complex systems consist- 
ing of up to hundreds of magnetic atoms within a sin- 
gle dot. [a, El 0] In this case, interparticle interactions 
along with anisotropy effects dominate the dynamics of 
the systems, and control the magnetization processes. [8| 
Moreover, since interdot exchange interactions are negli- 
gibly small, the dynamics of the ferromagnetic nanodot 
arrangements are strongly influenced by dipolar interdot 
interactions. 0, Hif 

Single molecule magnets (SMM's) consist of clusters 
of only a few magnetic ions, and are thus among the 
smallest and simplest nanomagnets, but are also well- 
characterized systems exhibiting magnetic hysteresis. 
In SMM's, the one-body tunnel picture of the magnetiza- 
tion mostly explains this phenomenon in the sense that 
the sequence of discrete steps in those curves provides 



evidence for resonant coherent quantum tunneling. [23, 
I29I l30| Nevertheless, this one-body tunnel model ne- 
glects intermolecular interactions, and is not always suf- 
ficient to explain the measured tunnel transitions. |3lll32l | 
A close examination of the magnetization curves re- 
veals fine structures which cannot be explained by that 
model. Wernsdorfer et al. suggest that these addi- 
tional steps are due to collective quantum processes, 
called spin-spin cross relaxation (SSCR), involving pairs 
of SMM's which are coupled by dipolar and/or exchange 
interactions. |3lL l32l| If dipolar and/or exchange interac- 
tions cooperate in the relaxation process, then one might 
hope to be able to better control such loss of magnetic 
memory. 

Analyzing the relaxation of the magnetization is dif- 
ficult for both SMM's and ferromagnetic nanodots. Be- 
sides dipolar interactions, many other factors may be in- 
volved in such processes. Geometric features, such as the 
shape and volume of the magnets, as well as the type of 
lattice on which they are placed, can directly influence 
the anisotropy barriers and the easy axis directions. In 
the case of SMM's, a quantum treatment has to be con- 
sidered to show that resonant tunneling of the magneti- 
zation results in the discrete steps appearing in the low 
temperature T magnetization curves. Although in many 
SMM's the intercluster exchange interactions are negligi- 
ble, as for ferromagnetic nanodots, in other SMM's, such 
interactions are comparable in strength to the dipolar 
interactions, Besides the quadratic Heisenberg and 
quadratic anisotropic intramolecular exchange interac- 
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tions, some SMM's are thought to contain intramolec- 
ular interactions of the Dzyaloshinskii-Moriya type. [3^1 
Additional higher order, anisotropic spin exchange inter- 
actions further complicate the problem. Therefore, by 
studying models that deal with each one of these factors 
separately, one hopes to simplify the problem, to build 
up gradually a more realistic system, and at the same 
time, to elucidate how each of these factors contributes 
to the magnetization process. 

With regard to SMM's, there have been recent 
approaches to the quantum dynamics of the low-T 
relaxation. [H |H IS HI El H3 Prokof'ev and 
Stamp assumed a single relaxation mode, [3] m which 
the dipolar and hypcrfinc fields are frozen unless an SMM 
flips its spin. Then by assuming the effective field around 
each SMM is that of randomly placed dipoles, they ob- 
tained an expression for the low-T decrease proportional 
to t p of the magnetization of each SMM from its fully 
magnetized state, 13 EJ ELwherep « 0.5 - 0.7, but 
p might be as large as 0.7.|34 l35t 13a [s^ . This proce- 
dure was restricted to very small deviations of the mag- 
netization from its saturated value, so it is not useful 
for studying the central portion of the hysteresis curves, 
for which the magnetization can be small. Moreover, as 
first argued for ferromagnets by Anderson, ^3] the spin- 
spin and spin-lattice relaxation times can be very differ- 
ent, so that such simple behavior is not expected. In 
fact, experiments on SMM's have shown that an expo- 
nential relaxation of the magnetization is consistent with 
the data. [33.13^1 so that as a minimum, one requires two 
distinct relaxation times for SMM's, which could be very 
different from one another. [4^ 

The most commonly studied model of spin dynam- 
ics containing two distinct relaxation parameters is the 
Landau-Lifshitz- Gilbert (LLG) equation. [43L liif Using 
the LLG equation, Kayali and Saslow (KS) investigated 
the hysteresis curves at T = for two-dimensional (2D) 
square arrays of 4 to 169 ferromagnetic nanodots subject 
to dipole-dipole interactions and a magnetic field applied 
in various directions within the array's xy plane. [45j They 
included anisotropy effects via an effective field propor- 
tional to the z-component of each dot's dipole moment. 
Earlier studies of square planar lattices of 9 to 36 fer- 
romagnetic dots were made by Stamps and Gamlev.jifij 
In addition, Zhang and Fredkin (ZF) studied the LLG 
model to obtain the zero-field time decay of the easy- 
axis magnetization of a three-dimensional (3D) cubic lat- 
tice ofl2xl2xl2 Stoner-Wohlfarth particles interacting 
with each other via dipole-dipole interactions. |l4j Since 
the size (or radius) of the Stoner-Wohlfarth particles was 
taken to be much less than the lattice constant, they 
could be treated as point-like magnetic moments, the 
classical analogue of SMM's. 

Here we study only the effects of the intermolecular 
dipole-dipole interactions upon the magnetization curves 
for an ensemble of iV c = 100 3D cubic crystals each 



containing N — 5 x 5 x 4 nanomagnets, all with the 
same magnetic moment. As in the ZF model of Stoner- 
Wohlfarth particles, we take the lattice parameter to be 
much greater than the nanomagnet size or radius. Except 
when a strong anisotropy field is present, we assume that 
there is no long-range order in the T regime of interest, 
so that in the absence of an external magnetic field, the 
magnetization of each nanomagnet crystal is essentially 
zero. We note that long-range ordering was claimed to 
exist in such a system with Ising spin anisotropy. [ItI lisj 
In our studies with a strong anisotropy field Ha, hys- 
teresis curves exhibiting a substantial zero-field magneti- 
zation were obtained for the applied magnetic induction 
B 1 1 Ha after the system had been fully magnetized by B. 
The strength of the dipole interactions is primarily deter- 
mined by the lattice constant, a, which we vary from 1.25 
nm to 2.5 nm. The dynamics of each nanomagnet are as- 
sumed to be given by the LLG equation, which includes 
precession and damping relaxation processes, the damp- 
ing coefficient a of which can also depend upon T. 
Then, the magnetic moment Mf of the ith. nanomagnet 
within the cth crystal of our ensemble obeys 
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where 7 = g[iB is the gyromagnetic ratio, M s = gfisS 
is the magnetic moment of an individual nanomagnet, 
and (B l c ) dip is the contribution to the effective magnetic 



induction Bf' cS at the ith nanomagnet within the cth 
crystal arising from dipole-dipole interactions between it 
and the other nanomagnets within the same crystal, 
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where the prime indicates that the j = i term is omit- 
ted. The second term of Eq. (1), the damping term, was 
first introduced by Landau and Lifshitz |43j | and later 
by Gilbert to give a phenomenological description of the 
relaxation of the magnetization. They did not derive it 
from first principles due to the enormous complexity of 
summarizing all of the relaxation processes into a single 
term. As noted above, when ferromagnetic interactions 
are present, a/7 is generally expected to be <C 1.01 
By extending to electronic spin systems the Wangsness- 
Bloch model of nuclear spin magnetic relaxation by mag- 
netic dipole coupling to a heat bath,^^| Fredkin and Ron 
showed that the damping term could be derived for large 
spin values and n = h"/H /ksT <C 1, where Ti and ks are 
Planck's constant divided by 27r and Boltzmanns' con- 
stant, respectively, and in our case H = _B^' eff .jH3| To 
the extent that electric quadrupole interactions could be 
neglected, a varies inversely with T for k <C 1, but de- 
pends upon k otherwise. |50| 
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More recently, a different derivation of the Gilbert 
damping term was derived from a spin Hamiltonian con- 
taining the interaction between the spin and the radiation 
field, which is induced by the precessing magnetization 
itself. 0, In that case, no explicit T dependence of 
a was given. We remark that rather complex explicit 
expressions for a for the different system of local spin 
moments arising from p — d kinetic-exchange coupling of 
the itinerant-spin subsystem in ferromagnetic semicon- 
ductor alloys have been given recently. |53l| In any event, 
the damping coefficient a at some T value must be de- 
termined experimentally for each system. 

In order to study the magnetization of ferromagnetic 
dots, KS used an extremely large value for the damping 
coefficient, a/7 = 0.6, a huge sweep rate, ~ 3000 
T/s, and a small maximum external induction B max = 
2fiQM s . In our studies of nanomagnet arrays, we used 
values of a/7 that varied from these values to values 12 
orders of magnitude smaller. Depending upon the a val- 
ues, we also varied the sweep rate AJ* from those values 
to the the much smaller ~ 10 -3 T/s, and varied B max 
from much larger values (2 T), comparable to those re- 
ported in SMM experiments. |23ll26| to those used by KS. 
Similarly, the lattice constants reported in the present 
work are in accordance with the near neighbor separation 
in the most extensively studied SMM crystals. Quantum 
processes within the individual SMM's will be treated in 
a separate publication. 

The present paper is organized follows. In Sec. II, we 
present our model system and a brief description of the 
overall calculation procedure that we followed. In Sec. 
Ill, we solve the LLG for each nanomagnet subject to 
both the external and the combined dipolar inductions. 
In Sec. IV, we present and discuss our main results for 
the magnetization curves, which are evaluated at vari- 
ous values of the sweep rate, T, a, the anisotropy field, 
and a/7. When spin anisotropy is present, we study the 
cases B \ \ Ha and B _L Ha- In Sec. V, we reproduce 
one of the KS magnetization curves for square 2D lat- 
tices, and vary some of their parameters to show that 
the results are almost independent of the sweep rate over 
the range ~ 300 — 6000 T/s. Analogously, we show that 
the anisotropy field does not affect significantly the mag- 
netization curves until its magnitude is comparable to 
Bmax/ /^o- By varying the lattice constant, we also show 
that the results of KS are very sensitive to the strength 
of the dipolar fields, which mainly govern the behavior of 
the magnetization of such systems. Finally, in Sec. VI, 
we summarize our main conclusions. 



MODEL SYSTEM 

In the present work we consider an ensemble of N c = 
100 cubic crystals (or configurations), each configuration 
consisting ofiV = 5x5x4 = 100 nanomagnets, each 



with ground state spin S = 5, which interact with one 
another only via dipolar interactions. Each of the N c sys- 
tem configurations c = 1 , . . . , N c is constructed to have 
a starting total magnetic moment M c w at B = 0. 
The hysteresis curves are obtained for each configuration, 
and these are then averaged over the iV c configurations. 
One then obtains the magnetization M. (B) curves, where 
M = (M c ) c /V is the configuration averaged magnetiza- 
tion, V is the crystal volume, and B = \B\. 

Ensemble of random spin configurations 

In order to proceed, we first find a large number N c 
of random spin configurations c of N = 100 spins, such 
that for each configuration, M c /M s « at B = and 
as T — > 00, where the total magnetic moment 

N=100 

M c (t,B) = M t&B). (4) 

i=l 

At the start of the iteration, we take t = 0, B = 0, 
and T — ► 00 in the absence of the dipole-dipole (or any 
other inter-spin) interactions for configuration c. Then 
we select those configurations for which \M C \/M S < 0.1, 
which we deem sufficiently close to M c sa 0. Our result- 
ing magnetization curves are based upon the average over 
N c = 100 configurations, each one containing N = 100 
similarly chosen nanomagnets. 

We reiterate that N is the number of nanomagnets in 
each configuration, and N c is the number of configura- 
tions studied. Although we have chosen both of these 
numbers to be 100 in order to obtain reliable statistics, 
N and iV c have completely different meanings. Finding 
many (N c ) configurations, each of which has an almost 
vanishing initial magnetization consumes a significant 
amount of computer time, especially if the number N of 
nanomagnets per configuration is not very large. How- 
ever, choosing a rather small number N of nanomagnets 
reduces the time required to calculate the dipolar field 
at each nanomagnet due to all of its neighbors, which 
must be calculated at each integration time step of the 
LLG equation, offsetting the large amount of computer 
time required to set up N c initially nearly-nonmagnetic 
configurations. 

Evolution of the magnetization versus field curves 

In this model one increases the external magnetic in- 
duction B in discrete steps AS, until B = B max , where 
B max = |.B max | has to be large enough to align ev- 
ery nanomagnet in its direction. How large B max has 
to be generally depends upon T, the field sweep rate 
= ^j-, the lattice parameter a, and the crystal 
struct urer|27j In addition, the steps \AB\ must be small 
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enough to give rise to numerically smooth A4(B) curves. 
We therefore set |AB| = B max /NB, where the number 
of steps Nb ^> 1 should be on the order of 10 3 . After 
each magnetic step, one allows each of the nanomagnets 
to relax for a fixed amount of time At, which is chosen to 
be sufficiently small that the nanomagnets do not reach 
equilibrium. Otherwise, in the absence of a sufficiently 
strong anisotropy field, no hysteresis would result. 

First, we choose one of our configurations c (e. g., 
c = 1) and set the moments of the nanomagnets equal 
to their values in this initially nonmagnetic configura- 
tion, {Mf =1 (t = 0,B = 0)}i=i,... <N . That is, just after 
we turn on the magnetic induction in the x direction by 
the amount B = AB, the nanomagnets have not yet 
precessed from their initial configuration. Then, we cal- 
culate the effective magnetic induction Bj~ 1,eS at each 
of the i = 1, . . . , N nanomagnets for c = 1. To do so, we 
must calculate the dipolar induction in Eq. (3) due to 
the presence of all the other nanomagnets. 

Then, we let each of the nanomagnets evolve in the 
presence of its effective magnetic induction for a chosen 
fixed time interval At. To do this accurately, we break 
At up into Nt intervals dt — At/N t . Obviously, this 
is extremely time consuming, because it is necessary to 
recalculate the effective induction at each nanomagnet 
after each time-integration step of width dt. Once the 
whole set of moments {M}(t = At, B = AB)}j=i jy is 
obtained, we proceed to calculate the configuration mag- 
netic moment, M c= i(At,AB), for this choice of fixed 
sweep rate, AB/At, from 



J2l 1 M?(t,B)eM-Wttt,B)} 



MJt, B) 

E-Liexp[-/3W?(t,B)] 
H?(t,B) = -B^ s {t,B).M?{t,B), 



(5) 
(6) 



by setting c = 1, t = At and B = AB, where B°' eS (t, B) 
is given from Eqs. J2J and 10, (3 = 1/kgT, and fc^ is 
Boltzmann's constant. Since (B 4 c ) dip as given by Eq. © 

in Bl' c (t, B) contains a self-fieldless single sum, there 
is no site overcounting in Eq. ©. 

We are interested in M C (B, ^f), a thermodynamic 
quantity that does not explicitly depend upon t. Al- 
though the Mf and B^' cS for each nanomagnet are eval- 
uated at each time step dt, the statistical weighting in 
Eqs. © and © is only evaluated at the end of each 
fixed interval At, which has a one-to-one correspondence 
with AB. Thus, this single-configuration average can 
be directly compared to those in different configurations 
after the same number of intervals. Moreover, since 
At = AB , M c (t, B) for our purpose can be writ- 

ten as M c ( sgTSt i B), which is effectively a function of 
BandAf. 

Next, we increase the external magnetic induction by 
another equal step AB, and let the nanomagnets precess 
during another equal time interval, At, under the action 



of the new effective induction. We continued increas- 
ing B in this equal step fashion a total of Nb times, 
until B = -B max - At this point, the incremental in- 
duction direction is reversed, setting B = B max — AB 
for the same time interval At, repeating the procedure 
2Nb times, until B = — B max . After that, we reverse 
the incremental induction direction once again, setting 
B = — -B m ax + AB for the same time interval At, etc., 
and continue Nb times until B = is reached, or un- 
til the configuration magnetization hysteretic loop (if it 
exists) is closed. Then, one repeats the entire procedure 
above described for each of the other N c — 1 configura- 
tions c = 2, . . . , N c , keeping the time intervals At and 
the subintervals dt constant for each step in each config- 
uration. Once all of the calculations for each of the N c 
configuration are finished, we average the results over the 
N c configurations, obtaining, 
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Then, the magnetization M. is easily calculated. Having 
tabulated M for every external magnetic induction step 
with fixed AJL ; T, a, M s , a, and -B ma x, we generate the 
magnetization curve A4(B) for this set of parameters. 



Variation of experimental parameters 

Unlike the parameters such as -Bmax and dt, which are 
details of the theoretical calculation, other parameters 
can in principle be varied in experiments in a variety of 
materials. Using the same initial dipole configurations 
we repeat the whole procedure with different values of a, 
AJ^ T, and a. The only parameters that can be experi- 
mentally varied in studies on a particular sample are A^ 
and T, since the other parameters are fixed. Neverthe- 
less, the possibility of setting the nanomagnets further 
apart by varying the composition of the non-magnetic 
ligand groups in SMM's, for example, justifies the study 
of the variation of a. Also, given that the damping term 
appearing in the LLG equation is phenomenological, and 
that in most cases a should be determined experimen- 
tally, we have also examined its variation. We note that 
a is expected to depend inversely upon T, unless T is suf- 
ficiently low that thermal processes no longer dominate 
the relaxation. We keep M s fixed. 



INTEGRATION OF THE LLG EQUATION FOR 
ONE NANOMAGNET 



The magnetic moment of each nanomagnet is obtained 
by numerically integrating the LLG equation. The time- 
evolution of one nanomagnet must be determined syn- 
chronously with all its neighbors in order to calculate 
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the dipolar induction acting on each of them at a given 
time. To solve the LLG equation for the ith nanomag- 
net in the cth crystal, we first rotate its coordinates at 
each time integration step such that B^' eS (t) || z(t). We 
then solve the resulting differential equations for either 
the coordinate spherical angles 9(t),4>(t), or the compo- 
nents of Mf(t), as shown in Appendix A. The quantity 
relevant to each spherical angle or component of Mf(t) 
is f* dr\B?' e (r)|, which explicitly involves the past his- 
tory of \B^' eS (t)\. In order to decrease the computation 
time, we approximate this integral for small time inte- 
gration steps dt = t — to -C to , 

fdr\B^ s {r)\ « \Bt> eS (t )\dt. (8) 

In order to assure numerical accuracy of our results 
for the greatly different experimental parameters stud- 
ied, we had to make appropriate choices for the numerical 
parameters used in the calculations, as discussed in the 
Appendix. Generally, calculations with slow sweep rates 
require correspondingly small a/7 values. For the calcu- 
lations leading to the results presented in Figs. 1, 2, and 
6-9, we take the numerical parameters dt = 1 x 10~ 4 s, 
B max = 2.0 T, N t = 1000, and N B = 500, 1000 and 4000, 
respectively, for the different sweep rates studied. For the 
calculations presented in Figs. 3-5, we take dt = 6x 10~ 12 
S) B max = 22.5 mT, N t = 1000, and N B = 1250. 

RESULTS AND DISCUSSION 

Effects of damping and maximum induction values 
on the hysteresis 

We first neglect any spin anisotropy effects. In Fig. 
1, we plotted the average over N c = 100 configurations 
of the normalized magnetization at the lattice constant 
a = 1.5 nm, sweep rate = 0.005 T/s, and tem- 
perature T — 0.7 K for the four weak damping rates 
a/7 = 3 x 10"", where n = 10, 11, 12, and 13. These re- 
sults appear respectively from left to right (right to left) 
in the upper (lower) part of Fig. 1. The magnetization 
curves show hysteresis for all four of these a values. For 
the smallest a value we studied, a/7 = 3.0 x 10~ 13 , the 
hysteresis only occurs for external induction magnitudes 
exceeding 3.0 T, observed by setting -B max above that 
value, which is well beyond the domain pictured in Fig. 
1. We also note that in Fig. 1, the central regions for 
\{M)/(NM S )\ < 0.8 are non-hysteretic. For each of these 
four parameter value choices, the initial curve describing 
the first increase of the average magnetization from essen- 
tially to its saturation value is indistinguishable from 
subsequent similar curves obtained after completing the 
full hysteresis paths. Hence, in this case, the main con- 
sequence of the choice of N c = 100 configurations is the 
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FIG. 1: Magnetization curves for 7V C = 100, a = 1.5 nm, 
%f = 0.005 From left to right for M > 0, a/7 = 3 x 10~ 10 
(dashed), 3 x 10" 11 (thin dark grey), 3 x 10" 12 (light grey), 
3 x 10~ 13 (thick black). The inset highlights the hysteretic 
region of the first three of these curves. 

improvement in the statistics, reducing the noise that re- 
mains most evident in the curves corresponding to the 
smallest a values. 

From the inset to Fig. 1, we see that although the 
height (in < M > /(NM S )) of the hysteretic region de- 
creases with decreasing a, the width (in B) of the hys- 
teretic region increases faster with decreasing a, so that 
the overall area of the hysteretic region increases with 
decreasing a. From a computational standpoint, for the 
parameter values studied in Fig. 1 , the smaller the value 
of a, the larger the required value of -Bmax to produce 
hysteresis. We also noticed that in these magnetization 
curves, the hysteresis sets in at the point of an abrupt 
change in slope in the initial curve, which describes the 
first increase of the average magnetization from to its 
saturation value. Moreover, we conclude that -B ma x must 
be chosen to guarantee that the system reaches satura- 
tion at B < -B max , because of the different nature of 
the hysteresis present in each curve. For example, in 
Fig. 1 the hysteresis can occur only after saturation, but 
with smaller a values, if the system has not saturated by 
-B = -B max , then the magnetization will keep increasing 
for a number of subsequent AB steps, even though the 
direction of AB (but not of B) has been reversed. 

Effect of temperature on the hysteresis 

Temperature-independent a 

We first investigate the role of temperature that arises 
only from the statistics, Eq. JSJ, and present our results 
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FIG. 2: Shown is the upper hysteretic region of the normalized 
magnetization curves at T — 0.7 K (grey) and T = 0.1 K 
(black, circles). The T- independent damping constants a/7 
are 3 x 10" 12 (a), 3 x 10 -11 (b), and 3 x 10" 10 (c). The other 
parameters are the same as in Fig. 1. 



for a T-independent a in Fig. 2. In this figure, we have 
replotted the inset of Fig. 1, excluding the curve for 
a/7 = 3 x 1CP 13 , for which the hysteresis occurred for 
B too large to display on the same plot. Otherwise, the 
parameters are the same as in Fig. 1, except that we 
have compared our results (grey curves) for T — 0.7 K 
shown in Fig. 1 with those (black curves and circles) for 
T = 0.1 K. Since the evolution of the magnetization with 
B in this model is independent of T, we note from Fig. 
2 that the departures of the magnetization curves from 
the points of saturation are the same at both T values, so 
that the widths (in B) of the hysteretic regions are nearly 
the same. However, the height in < M > /(NM S ) of each 
hysteretic region decreases strongly with decreasing T, so 
that the overall area of each hysteretic region decreases 
with decreasing T. This particular result is in strong 
contrast to the existing experimental results on SMM's. 
Nevertheless, our results are reasonable from the point 
of view of the LLG equation and the way T enters our 
calculation. We reiterate that we have so far neglected 
quantum and spin anisotropy effects, the latter of which 
will be discussed in the following. 

We remark that in Fig. 2, T only enters into the equa- 
tions of motion when the average magnetic moment is 
evaluated from Eq. JSJ. As for the Brillouin function 
that describes the magnetization of a paramagnet in the 
absence of the dipolc interactions, the initial slope of the 
magnetization at low B increases as T is lowered. This 
increases the alignment of the moment of each nanomag- 
net at decreasing T, so that the dipole-dipole interactions 
tend to be maximized, enhancing the effect. 



Temperature-dependent a 

We now consider the effect of the temperature depen- 
dence of the damping constant a upon the magnetic hys- 
teresis, focussing upon the case of correspondingly fixed 
very high sweep and damping rates. We assume that 
our choice of spin value, S = 5 for each nanomagnet, 
satisfies S ^s> 1. In this limit, Fredkin and Ron showed 
that the damping of the nuclear spin precession by mag- 
netic dipole coupling to a heat bath, as derived under the 
assumption of spin-orbit factorization by Wangsness and 
Bloch, could be readily extended to the spins in magnetic 
systems. H<J|§3 For S > 1, they found 

<*{T)h « % (9) 

where T Q = 2h<b\ l (1 - e~ K )S 2 /k B K,]j$ and $ 1X is a rate 
constant (with units of s _1 ), the expression for which is a 
complicated orbital integral arising from the interaction 
of the local spin with its surrounding molecular electronic 
orbitals in second-order perturbation theory. |4flj and n — 
h~fB cff /(k B T). For k < 1, T -> 2h<t>{ 1 S 2 /k B , which can 
be taken to be independent of T and B cS , so that a oc 
T _1 , but for k 1, a oc l/B cS , which would completely 
change its effect. Here we only consider the case k <C 1, 
for which Eq. @ holds for constant To. We note that, 
as in Figs. 1 and 2, T also affects the results for the 
magnetization from the statistics, Eq. JSJ. 

In Fig. 3, we have shown our results, averaged over 
N c = 100 configurations, of the normalized magnetiza- 
tion as a function of B in mT, for a = 1.5 nm, = 3000 
T/s, a(T)/ 7 = Tq/T, T = 0.3K, and T — 5 K. For 
the calculations presented in this figure, we used the nu- 
merical parameters dt = 6 x 10 -12 s, B max — 22.5 mT, 
N t = 1000, and N B = 1250. Note that although a has the 
same value as in Figs. 1 and 2, the sweep and damping 
(a(T)/j = 0.06) rates are six and at least eight orders of 
magnitude larger than in those figures. For these param- 
eters, there are three regions of hysteresis in the pictured 
magnetization curve. The left inset is an enlargement of 
the upper hysteretic region, the mirror image of which 
occurs in the lower region of the pictured magnetization 
curve. In contrast to the behavior shown in Figs. 1 and 
2, at the top of the upper hysteretic region, the mag- 
netization does not rise abruptly to its saturation value, 
but first goes through an extended non-hysteretic region. 
In addition, there is a central hysteretic region, an en- 
largement of which is shown in the right inset, along 
with an enlargement of the same central hysteretic re- 
gion obtained at T = 0.25 K with the same set of pa- 
rameters. We note that at T — 5 K, the onset magneti- 
zation averaged over N c = 100 configurations, pictured 
by the thin curve in the lower portion of the right inset, 
does not coincide with the thick curve corresponding to 
the central hysteresis loop region obtained subsequently 
to the attainment of the saturation value by the mag- 
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FIG. 3: The magnetization curves for N c = 100 at T = 5 K, 
a = 1.5 nm, and %f = 3000 T/s with B max = 22.5 mT and 
«( T )/7 = T o/T for h gf i B B efi ' /(k B T) < 1 and T = 0.3 K 
are shown. |50B Left inset: details of the upper portion of the 
curve. Right inset: details of the central hysteretic portion 
of the curve shown, along with the central portion of the cor- 
responding curve at T — 0.25 K. The thin curves beginning 
near to the origin represent the magnetization onsets. These 
curves are offset for clarity, with the scales on the right (left) 
sides corresponding to the lower (upper) curves, respectively. 




-1.0 -0.5 0.0 0.5 1.0 
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FIG. 4: The central loop and starting magnetization curves 
for two separate configurations, each with N c = 1 (open and 
filled circles) at T = 10 K, a = 1.5 nm, and %f = 3000 T/s 
with a{T)H = To/T for hg^ B B cf! /k B T < 1 and T = 0.3 K 
are shown. |5Cl The thin grey and thick black curves represent 
the identical behaviors of the central hysteretic loop portion of 
the magnetization for the same two configurations obtained 
after saturation. The arrows indicate the direction of the 
magnetization hysteresis. Here B max = 22.5 mT. See text. 



netization. In addition, we note that the thick central 
hysteresis loop exhibits reproducible oscillations with B- 
independent frequency / at T = 5K, which oscillations 
have disappeared at the lower T = 0.25 K value, for 
which a(T)/"f = 1.2, pictured in the upper portion of 
the right inset. 

In order to investigate further the differences between 
the starting magnetization curves and the curves ob- 
tained subsequent to saturation, in Fig. 4, we have 
shown the corresponding central hysteresis loop portion 
of the magnetization obtained for two individual config- 
urations, using the same experimental and numerical pa- 
rameters as in Fig. 3, except that T = 10 K, for which 
a(T)/7 = 0.03. As in Fig. 3, T enters both through 
the statistical averaging and through the damping, a(T). 
In Fig. 4, the solid and open circles correspond to the 
starting magnetizations of the two configurations, and 
the coincident thick black and thin light grey curves cor- 
respond to the central hysteresis loop region of their re- 
spective magnetization curves obtained after saturation. 
Note that after the initial noisy regions, the starting mag- 
netizations for these two configurations exhibit compa- 
rably large amplitude oscillations at the frequency f/2, 
the phases of which are very different. However, after the 
attainment of the saturation magnetization, these large 
amplitude oscillations are absent, and replaced by smaller 
amplitude oscillations at the frequency /, which are sim- 
ilar to the oscillations present in our results obtained at 



T = 5 K shown in the lower curves in the right inset of 
Fig. 3. We note that in the first oscillation present on 
both sides of the central post-saturation hysteresis loops 
obtained with these parameters at T = 5 and 10 K show 
additional small amplitude, higher frequency oscillations, 
which may be higher harmonics of /. In addition, the 
amplitudes of the fifth and sixth oscillations are larger at 
10 K in Fig. 4 than at 5K in the lower right inset of Fig. 
3. 

We remark that the large amplitude oscillations 
present in the starting magnetizations shown in Fig. 4 
are absent in Fig. 3. This occurs due to the randomness 
of the oscillation phases, which is averaged out in the 
N c = 100 configurations studied in Fig. 2. 

From Fig. 4, we therefore conclude that our starting 
configurations that were chosen to have \M\/M S < 0.1, 
appropriate for SMM's, lead to starting magnetization 
curves that are very different from those that start at 
the saturation magnetization, but are subsequently iden- 
tical. That is, after the attainment of saturation, all 
configurations arc identical. 



External field directed towards the crystal corners with a(T) 

We now consider the 3D case of the external magnetic 
induction directed from the crystal center to one of its 
corners, B = B(x+y)/y/2, the (110) direction. In Fig. 5, 
we show the resulting central hysteresis region obtained 
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FIG. 5: The central loops (solid curves) of the magnetization 
curve for N c = 50, TV = 5 x 5 x 4 = 100, with B along the 
[110] direction [B||(x+yW2] at T = 10 K, with a = 1.5 nm, 
B max = 22.5 mT, and %f = 3000 T/s with a(T)h = T /T, 
and To = 0.3 K. The dashed curve is the starting magnetiza- 
tion curve. The arrows indicate the direction of the hysteresis. 
Inset: the full magnetization curve. See text. 

from our calculations for N c = 50, N — 5 x 5 x 4, T = 10 
K, a = 1.5 nm, and %f = 3000 T/s with a(T)/~f 
0.03. In this case, it is sufficient to set B max = 22.5 
mT, which leads to full saturation. We note that for this 
field direction, a small (-6 mT< B <6 mT) hysteresis 
region appears on either side of the origin, which is rather 
central to the full magnetization curve, but vanishes over 
a small region close to the origin. There are also tiny 
hysteresis regions near to saturation that appear as dots 
in the inset depicting the full magnetization curve. 

The nearly central hysteretic regions shown in Fig. 5 
exhibit reproducible jumps at specific B values, similar 
to those observed at low T in SMM's. However, we note 
that in this figure, we have taken T = 10 K, and have 
used a classical spin model. We also note that we have 
used a rather small sample (N = 100) with a fast sweep 
rate and a large damping coefficient in our calculations, 
and caution that such behavior might not be present in 
large single crystals, especially with much slower sweep 
rates. Nevertheless, this figure demonstrates that steps 
in the magnetization do not necessarily have a quantum 
origin, and that the sample shape can lead to unusual 
hysteresis effects. 



Effect of sweep rate on the hysteresis 

From the curves obtained using the same numerical 
parameters as in Fig. 1 for different induction sweep 
rates at a fixed, small damping rate in Fig. 6, it is clear 




B(T) 



FIG. 6: Hysteretic region of M(B) at 0.7 K, a/7 = 3x 10" 12 , 
and a = 1.5 nm, for the sweep rates = 0.04 T/s (thin 
black), 0.02 T/s (dark grey), and 0.005 T/s (thick light grey). 
The inset shows the entire curves. 

that stronger hysteresis is found for higher sweep rates, 
in agreement with experiments on a variety of nanomag- 
nets, including SMM's. This shows that the reversibility 
of the process depends on how close to equilibrium the 
sweep rate allows the nanomagnet spins to reach. That is, 
although for different sweep rates the external induction 
is increased by the same amount AB, at higher sweep 
rates, the time At allowed for the nanomagnets to evolve 
towards equilibrium is less. This makes the process less 
reversible and the hysteresis loops larger. 

We also note that at the much higher sweep and damp- 
ing rates studied in Figs. 3, 4, the magnetization also 
exhibits a central hysteretic region, which exhibits oscil- 
lations at T values not too low and/or damping constants 
not too large. 

Effect of lattice constant on the hysteresis 

In Fig. 7, we show hysteresis curves for two different 
values of the lattice constant a, obtained using the 
same numerical parameters as in Fig. 1. For weaker 
dipole-dipole interactions (larger a), the rise in the 
magnetization is steeper with increasing B, and the 
rapid decrease in the magnetization from its saturation 
value upon decreasing B occurs at a smaller value of 
\B\. Furthermore, we shall see that dipolar interactions 
do not promote hysteresis in these systems, but suppress 
it. Actually, the same conclusion was found recently 
for magnetic nanoparticles in the framework of the 
generalized mean-field approximation. [55j 

This peculiar hysteresis is easily understood by ana- 
lyzing the LLG equation. If the nanomagnet magnetiza- 
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region of the Ai(B) curves. We remind the reader of our 
intent to focus upon the effects of the dipole-dipole inter- 
actions, whereas the most important features of SMM's 
involved in their low-T relaxation of the magnetization 
are generally thought to be their quantum structure and 
magnetic anisotropy. Nevertheless, for this entirely clas- 
sical and magnetically isotropic system, we are indeed 
finding hysteretic curves. In addition, the sweep rates in 
Figs. 1, 2, 6, and 7 are comparable to those used in ex- 
perimental SMM studies. At much larger sweep rates, 
such as were studied in Figs. 3-5, an hysteretic cen- 
tral region was found. However, the sizes and T depen- 
dencies of these hysteretic regions were still respectively 
much smaller and qualitatively different than observed in 
SMM's. 



FIG. 7: Magnetization curves for lattice constants a — 1.5 
nm (grey) and a = 2.5 nm (black). %f = 0.04 T/s, a = 
3 x 10~ 12 7, T = 0.7 K. 

tion Mf is parallel to its local magnetic induction Bf ,c , 
dMf jdt = 0, as it will remain thereafter, so that Mf has 
reached equilibrium. The only chance the system has 
to decrease its magnetization from its saturation value 
is through the combined weak dipolar induction, which 
strengthens with decreasing lattice parameter a. The 
dipolar induction can oppose the system from remain- 
ing completely magnetized, since it has small, but non- 
vanishing components. Therefore, even when B reaches 
its maximum (finite) amplitude B max and the misalign- 
ments of each Mf with B are negligible, dynamic equi- 
librium will not generally have been attained due to the 
limited time allowed for relaxation before the next change 
in B. There will remain a slight deviation between the 
directions of the Bf ,cS and the Mf due to the presence 
of the Bf ,dvp \ which is especially important when B de- 
creases from -B max - 

Of course, it is harder to decrease M. at the very be- 
ginning of the induction cycle. This is precisely the cause 
of the hysteretic behavior, given that changes in Mf are 
proportional to \Mf x B^' eS | , which nearly vanishes when 
the direction of the incremental induction has just been 
reversed. We conclude then that the smaller the lat- 
tice parameter (the stronger the dipolar induction), the 
greater the deviation of Mf from the direction of B- ff . 
Hence, the easier it is to decrease M, making the magne- 
tization curve less hysteretic. This is shown in Fig. 7, in 
which the magnetization curves resemble those obtained 
for Mn4 SMM's. Those data show an abrupt decrease 
in M. at nearly zero external induction that is not evident 
in the magnetization curves of other SMM's. 

It is important to note that the curves in Figs. 1-7 do 
not show the strong hysteresis observed experimentally 
in most SMM's, which is especially large in the central 



Effect of spin anisotropy upon the hysteresis 

It is straightforward to generalize our model to include 
some of the effects of magnetic spin anisotropy. Here we 
assume the nanomagnets contain sufficiently many spins 
that their quantum nature can be neglected. We note 
that SMM's at low T values behave as quantum entities, 
because of the small number of spins in each nanomag- 
net. In those systems, most workers have assumed the 
in addition to the isotropic Heisenberg and Zeeman in- 
teractions, the magnetic anisotropy terms could also be 
written in terms of components of the global spin op- 
erator S, with the overall dominant terms often writ- 
ten as —DSl — E(S% — S^-Hil However, portions suffi- 
ciently large for model comparison of the low-T magne- 
tization curves of two Fe? SMM dimers have been stud- 
ied experimentally. [57l l58| In neither antiferromagnetic 
dimer case was any evidence for either of those types of 
spin anisotropy present. 59] In contrast, in one of those 
cases, strong evidence for a substantial amount of local, 
single-ion spin anisotropy, in which the individual spins 
within a dimer align relative to the dimer axis, is present 
in the data. [H?], In addition, the global anisotr opy 
in the ferromagnetic SMM Mng is extremely weak. |60j 
Since the precise quantum nature of more complicated 
SMM's appears therefore to be poorly understood, we 
shall investigate the quantum features of the magnetic 
hysteresis curves in SMM's, including some effects of lo- 
cal spin anisotropy, in a subsequent publication. |54j 

We therefore restrict our investigations of the role of 
magnetic anisotropy upon the magnetization curves of 
arrays of nanomagnets to the simplest classical model of 
spin anisotropy, 

B o,en = B + B C d . p + ^ QHAi (10) 

where we take B = Bx and studied the cases Ha = Hax 
and Ha = H az. This is the 3D analogue of the model 
studied by KS. 45J In this model, the magnetic anisotropy 



10 




B(T) 



FIG. 8: Parallel 3D magnetization curves including differ- 
ent anisotropy field Ha = Hax strengths, with the exter- 
nal induction B\\Ha- HoHa — (thin black), 0.2 T (dark 
grey), and 1.0 T (thick dashed), respectively. For each curve, 
iV c = 100, Af = 5x5x4 = 100, a/j = 3 x 10~ 12 , a = 1.5 
nm, %f = 0.04 T/s, T = 0.7 K, and B max = 2.0 T. 



of each of the nanomagnets points in the same direction, 
and in our finite sized crystal consisting of 5 x 5 x 4 nano- 
magnets on a cubic lattice, our chosen direction is one 
of the most general ones. We first performed two stud- 
ies of the magnetic hysteresis in this model, for which 
the anisotropy field Ha is directed respectively along 
(100), H-B, and (001), _L B, and our results are shown 
in Figs. 8-9, respectively. For both anisotropy field di- 
rections, we take N c = 100, JV = 5 x 5 x 4 = 100, 
a/ 7 = 3 x 10~ 12 , a = 1.5 nm, = 0.04 T/s, T = 0.7 
K, and B max = 2.0 T. The sweep rates used in Figs. 8-9 
are slightly faster than those used in SMM experiments 
but much slower than those used in the calculations of 
KS. Since a = 1.5 nm in these curves, these curves also 
represent the strongest realistic dipolar interaction we 
studied. 

In Fig. 8, we show the portions of the parallel magne- 
tization curves with B||ff^||x, that exhibit the resulting 
regions of magnetic hysteresis for three Ha values. For 
the Ho Ha = 0,0.2, and 1.0 T values shown, all three 
curves are hysteretic, but the two lower Ha values do 
not lead to a central hysteresis region. Nevertheless, the 
largest anisotropy value, Ha = 1.0 T, leads to a strong 
central hysteresis. We remark that the trends shown in 
Fig. 7 are rather different from those obtained for a single 
magnetic particle with magnetic anisotropy. 

In Fig. 9, we show the portions of the 3D perpen- 
dicular magnetization curves exhibiting the resulting re- 
gions of magnetic hysteresis for the five anisotropy fields 
(i H A = 0, 1 mT, 12 mT, 0.5 T, and 1.0 T, with the 
magnetic induction B\\x _L iJ^Hz. In each case, hys- 




FIG. 9: Upper region of the 3D perpendicular magnetization 
curves with the external induction B — Bx _L Ha = Haz, 
for different values of Ha- Curves (a)-(e) correspond to 
HoHa = 0, 1 X 10 -3 , 1.2 x 10" 2 , 0.5, 1.0 T, respectively. 
The other parameters are the same as in Fig. 7. The arrows 
indicate the directions of the field sweeps. 



teresis occurs near to magnetic saturation, but is absent 
in the central region for small magnetic induction. At 
HqHa = 1.0 T, this is distinctly different from the large 
central hysteretic region observed for parallel anisotropy. 
Note that the slope AM / dB at small B is non-monotonic 
with increasing Ha, as it has a minimum at curve (c), 
corresponding to /j.qHa — 12 mT. 

Thus, we conclude that it is possible to obtain a cen- 
tral hysteresis region using this classical model of dipo- 
lar interactions with constant spin anisotropy. How- 
ever, our results suggest that such central hysteresis re- 
gions only arise for the magnetic induction parallel to 
the spin anisotropy direction, and for sufficiently strong 
anisotropy fields, H A > H^ n , where 1.0 T> /i -ff™ in > 
0.2 T. 



DIPOLAR INTERACTION, INDUCTION SWEEP 
RATE, AND ANISOTROPY DEPENDENCIES 
FOR A 2D SYSTEM 

To estimate the importance of the dipolar induction 
(especially when it becomes comparable to the external 
induction), the anisotropy and the sweep rate, we have 
reproduced one of the 2D calculations of KS.|45j The KS 
calculation we chose to reproduce was pictured in their 
Fig. 2(i), and is shown here as the left panel of Fig. 10. 
Then, we changed some experimental parameters to see 
how the results depend on the anisotropy strength, sweep 
rate, and lattice parameter. 

Our calculations for a cubic lattice consisting of four 
25-particle layers differ from those of KS in many 
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ways.[45j They used a 2D square lattices of cylindrical 
nanodots (here, we take their 5x5 lattice with external 
induction aligned along an array's diagonal), included 
a shape-dependent anisotropy field perpendicular to the 
lattice, performed their calculations at T = 0, used a 
much larger damping constant than we generally did it 
for 3D systems, and did not average their results over 
an ensemble of 2D samples, because such systems do not 
show variations in the resulting hysteresis loops for dif- 
ferent initial states. Nevertheless, we both integrated the 
LLG equation using the Runge-Kutta algorithm, and sur- 
prisingly, KS's system turned out to be very sensitive to 
the dipolar field strength. The effective induction they 
considered can be written as 

B c, cS = B + B e dip + fM)HA ^ (n) 

For lattice constant a — 1.5 nm, spin S = 5, and 
V/a 3 — 0.5, where V is the volume of the nanomag- 
net, the saturation magnetization is M s ps 55 Oe. Then, 
they took the dimensionless dt = 5 x 1CP 3 , which im- 
plies a real time interval dt = 5.17 x 10~ 12 s. If the 
system evolves during 700 time steps dt for some fixed 
value of B, then B is changed every At w 3.62 x 10~ 9 
s. On the other hand, KS chose a maximum exter- 
nal induction £? max = 2(j,qM s ps 1.1 x 10 -2 T. In ad- 
dition, they took fixed induction steps of magnitude 
AB = 2 x 10- 3 noM a fa 1.1 x 10~ 5 T. Therefore, we 
estimate their resulting sweep rate to be yM- « 3 x 10 3 
T/s, as in our 3D results shown in Figs. 3, 4. 

In the absence of any specific information, we then had 
to induce the value of the anisotropy field that KS used 
to obtain their figure. Fortunately, as discussed in the 
following, the results are rather insensitive to it, unless 
Ha becomes comparable to -B m ax/Mo- in the right panel 
of Fig. 10, our 2D calculation with [IqHa — 0-75 mT are 
shown, and by comparing that figure with Fig. 2(i) of 
KS pictured in the left panel of Fig. 10, we see that the 
agreement is remarkably good. 

Very recently, Takagaki and Ploog (TP) used a fourth- 
order Runge-Kutta procedure to integrate the LLG equa- 
tions with N x N 2D nanomagnet lattices with magnetic 
anisotropy and dipole-dipole interactions . |6 1) They used 
a fixed time interval dt = O.lh/ (jM s ), 20 times as fast 
as that used by KS,0| and continued interating until 
no further changes in the nanomagnet spin configura- 
tions were obtained. They obtained results for TV = 5 
which they described as considerably different from those 
of KS, with a somewhat different magnetization loop and 
a larger area of the hysteretic regions. |61| Although they 
claimed that their fourth-order procedure was intrinsi- 
cally more accurate than the second-order one used by 
KS and by us, the fact that we obtained the excellent 
agreement pictured in Fig. 10 with one of the N = 5 
results of KS suggests that the procedure used by TP 
might have been less accurate than they claimed. |61| 
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FIG. 10: (left) Hysteresis loop M(B) in units of M s , for a 
weakly coupled array of 5 x 5 ferromagnetic nanodots in a 
square lattice on the xy plane, from Fig. 2(i) of KS. The 
external induction is applied along the array diagonal (45° 
from the x axis) . 45] . (right) Our results calculated for N c = 1 
with 5x5 nanomagnets on a square lattice, a/y = 0.6, T = 
OK, = 3000 T/s, h Ha = 7.5 x 10~ 4 T, B = B(x+y)/V2. 
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FIG. 11: Hysteresis loops for different strengths of Ha for 
5x5 nanomagnets on a square lattice with N c = 1. S = 5, 
T = K, a = 1.5 nm, = 3000 T/s, a/7 = 0.6, B = B{x + 
y)/v2. The thin grey and thick black curves with no Ha = 
0,0.75 mT, respectively, are nearly indistinguishable. The 
small grey circles and dashed curves correspond to [ioHa = 
4.0, 5.5 mT, respectively. The inset shows the entire curves, 
which are symmetric with respect to the origin. 



Anisotropy field dependence of the hysteresis 

We first investigated the effects of changing the 
strength of the anisotropy fields, and presented our re- 
sults in Fig. 11. The most important issue about the 
results shown in Fig. 11 is the fact that the curve ob- 
tained by KS (the left panel of Fig. 10) is basically inde- 
pendent of the anisotropy field Ha for sufficiently small 
Ha- That is, there are no essential differences between 
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that curve reproduced in the right panel of Fig. 10 with 
Ho Ha = 7.5 mT, and the one with Ha — 0. Strong 
deviations from these essentially identical curves appear 
for /j.qHa > 4 mT, however. Since identical behavior 
is obtained without any anisotropy, this implies that all 
hysteretic features (including the stepped magnetization 
and demagnetization) are due to the dipolar interaction. 
Ha becomes important only when it is comparable to 
B max / Ho and tends to close the hysteresis loops, starting 
from the lower and upper loops. 

We note that by comparing Fig. 11 with Fig. 9, the de- 
tails of the hysteresis obtained with Ha = for B along 
the (110) direction are different in 3D and 2D samples. 
The hysteresis is much larger in the 2D case pictured in 
Fig. 11, and has a large loop in the central region that 
does not vanish at the origin, plus large loops that ex- 
tend up to saturation. In the 3D case constructed from 
four 2D planes each equivalent to that used in the calcu- 
lation shown in Fig. 11, the magnitude of the hysteresis 
is reduced and its details have been greatly altered. 



Induction sweep rate dependence of the hysteresis 
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FIG. 12: Hysteresis loops for different induction sweep rates 
with 5x5 nanomagnets on a square lattice with N c = 1. 
S = 5, T = K, a = 1.5 nm, fi H A = 0.75 mT, a/7 = 
0.6. The dashed grey, thick black, and light solid grey curves 
correspond to = 300, 1500, 6000 T/s, respectively. The 
inset shows the entire curves. B = B(x + y)/\/2. 



In Fig. 12, we show our results for a single configu- 
ration of a square 2D lattice with JV = 5 x 5 for differ- 
ent sweep rates, keeping the other parameters fixed at 
HoH A = 0.75 mT, a/7 = 0.6, a = 1.5 nm, S = 5, T = 0, 
and B = B(x + y)j\f2. From Fig. 12, we note that 
the hysteresis is nearly independent of induction sweep 
rate over the range 300 to 6000 T/s, distinctly different 
from the strong dependence found in 3D systems shown 
in Fig. 5. 



Lattice parameter dependence of the hysteresis 

In Fig. 13, we have illustrated the effect of the lattice 
constant a upon the hysteresis. In this figure, we kept the 
other parameters fixed at S = 5, T = 0, %f = 3000 T/s, 
HoH A = 0.75 mT, a/7 = 0.6, and B = B(x + y)/V2. 
As a is varied from 2.0 to 1.25 nm, the upper portions 
of the hysteresis curves appear from left to right, respec- 
tively. From Fig. 13, it is readily seen that the magne- 
tization curves are very sensitive to a and hence to the 
strength of the dipolar interaction, which is proportional 
to a~ 3 . Our results for a=2.5 nm exhibit a smaller hys- 
teresis shifted further to the left, and all indications of 
steps have disappeared. Although not shown in Fig. 13, 
as a is increased further to 3.0 nm, the hysteresis almost 
disappears entirely. We deduce that stronger dipolar in- 
teractions (smaller a) result in larger hysteresis loops con- 
taining increased widths of additional steps. 

We then infer that contrary to the conclusion found 
for the 3D systems (based upon much smaller damping 
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FIG. 13: Hysteresis loops for lattice parameters a — 2.5 nm 
(solid black), a — 2.0 nm (dashed black), a = 1.5 nm (solid 
grey), and a — 1.25 nm (dot-dashed black), for 5x5 nano- 
magnets on a square lattice with N c = 1. S = 5, T — K, 
%f = 3000 T/s, (i H A = 7.5 x 10~ 4 T, a/7 = 0.6. The inset 
shows the entire curves. B = B(x + y)/V2. 



coefficients and much slower sweep rates), the dipolar in- 
teractions promote a hysteretic behavior in this 2D sys- 
tem. 
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SUMMARY AND CONCLUSIONS 

We first found N c = 100 sample configurations with 
an overall magnetization close to 0. We then solved the 
Landau-Lifshitz-Gilbert equation for a 3D cubic lattice of 
N = 5x5x4 nanomagnets, subject to dipole-dipole inter- 
actions and spin anisotropy. These results should be rel- 
evant for an array of Stoner-Wolfarth nanomagnets, and 
to some extent, single molecule magnets, although the 
quantum nature of the latter has so far been neglected. 
In the absence of spin anisotropy, we varied the magnetic 
induction sweep rate ^j, the damping constant a, the 
lattice constant a, and the temperature T. We also con- 
sidered the effects of a T-dependent damping constant 
of the form a(T)/7 = Tq/T suggested by Fredkin and 
Ron. For slow sweep rates and small a relevant for ex- 
perimental studies on single molecule magnets, magnetic 
hysteresis appears in the regions of the magnetization 
curves near to saturation, the area and onset magnetic 
induction strength of which increases with decreasing a 
and increasing sweep rate. With decreasing T, the onset 
magnetization magnitude of the hysteretic regions near 
to saturation decreases. With decreasing a correspond- 
ing to increased dipole-dipole interaction strengths, the 
onset of the hysteresis regions near to saturation appears 
at increasing magnetic induction magnitude. 

At much larger sweep rates and damping constants, the 
magnetization curves attain saturation at much smaller 
applied magnetic induction strengths. The hysteretic re- 
gions just below saturation have moved somewhat be- 
low saturation, and a new central hysteretic region ap- 
pears. As one follows the magnetization curve for a sin- 
gle configuration, the starting curve exhibits oscillations 
at a rather constant (magnetic induction independent) 
frequency //2, but the phase of the magnetization oscil- 
lations is a random function of the configuration. After 
the attainment of magnetic saturation, this central hys- 
teretic region exhibits oscillations at /, twice that fre- 
quency, possibly with weak higher harmonics, for T not 
too low, which are independent of the configuration. 

When the applied magnetic induction is in the (110) 
direction (from the sample center to one of its corners), 
magnetic hysteresis exhibiting steps and jumps appears 
within the central region, but vanishes at and very near 
to the origin. Although these step-like features are sug- 
gestive of the behavior seen in single molecule magnets, 
they are present at rather high T values, as they arise 
from the classical sample shape effects. 

In the presence of the magnetic anisotropy field Ha, 
an applied magnetic induction parallel to the anisotropy 
axis leads to a large central hysteresis region, provided 
that the magnitude of the spin anisotropy is sufficiently 
large. For the applied magnetic induction perpendicular 
to the magnetic anisotropy, no central hysteresis region 
is present, although a small amount of hysteresis near to 



saturation persists for sufficiently small spin anisotropy 
strength, and the slope of the magnetization curve at the 
origin is non-monotonic, exhibiting a maximum at a par- 
ticular small value of the spin anisotropy strength. These 
effects for the spin anisotropy are qualitatively in agree- 
ment with those in many types of arrays of nanomagnets, 
including single molecule magnets. 

As a test of our numerical procedure, we studied the 
simplified 5x5 2D square lattice with a perpendicular 
spin anisotropy field Ha using the same procedure, and 
for a particular set of parameters, obtained quantitative 
agreement with a hysteresis curve obtained for that sys- 
tem by Kayali and Saslow. We showed that their hys- 
teresis curve is basically independent of Ha until ^qHa 
is on the order of the external induction. We also demon- 
strated that the magnetic hysteresis does not depend sig- 
nificantly upon the magnetic induction sweep rate, as op- 
posed to the dependence we found in our 3D system. In 
addition, we found that the magnetization of the 2D sys- 
tem is very sensitive to variations in the lattice parameter 
a. Finally, we noticed that although dipolar interactions 
also oppose the magnetization process in 2D systems, in- 
creasing the onset magnetic induction strength for the 
attainment of saturation as in 3D systems, they increase 
the area of the hysteresis, a behavior opposite to that 
found for the 3D system with a much smaller damping 
coefficient and much slower sweep rate. 

We expect our results to be relevant to the magneti- 
zation processes in a variety of nanomagnet arrays, es- 
pecially those approximating arrays of Stoner-Wolfarth 
particles. In addition, some of the features we obtained 
should be relevant to single molecule magnets, although 
the temperature dependence of the effects is not in agree- 
ment with experiments on those materials. Further stud- 
ies of the magnetic hysteresis using quantum models of 
the nanomagnets and their various anisotropy types is 
warranted, and will be addressed subsequently. |5J| 
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APPENDIX A 

We rotate our reference frame at every integration time 
step in such a way that Z??' eff is along the z axis. In this 
case, we can easily solve the LLG equation, Eq. (1). 
For simplicity of notation, we drop the subscripts i and 
superscripts c, and remember that we are describing the 
precession of the ith nanomagnet in the cth crystal. We 
define the axes to describe the magnetization direction 
of this particular nanomagnet, M, 8, and 0, where <fi = 



14 



M x 0, and then write 

B eff 



= B z z = B z (Mcos6 - 0smt 
= MB M + 6B g . 



(12) 



Since the magnitude of the dipole moment M s is con- 
served, in spherical coordinates Eq.(l) leads to 



dM ~d6 2 . M 

— — = 0— +4>sm0-i- 
dt dt dt 

= 8aBg + 4>jBg. 



Finally, from 



d9 
lit 



4>(to + dt) 
9 {t + dt) 



(13) 

(14) 
(15) 

(16) 

0(to)-a\B cS (to)\sm[9(t o )}dt. (17) 



= -a|B cff |sin( 

ffi 



ft = " 
we obtain for a very small time interval dt, 



<Xi )-7|B Gff (W* 



These equations were used in our numerical calculations. 
In order to relate the angles to measurable quantities, 
however, we note that it is possible to integrate Eqs. i|14|) 
and l|15f) exactly, obtaining 



cos 



tanhltanh 1 (cos[0(io)] 



6(t) 



+a / dr\B cn (r)\ 

J t 

</>(t) = 0(to) -if dr\B ef{ ( 



(18) 
(19) 



which is equivalent to that obtained using a somewhat 
different technique. We note that by expanding Eqs. 
(fT5f) and [EH to leading order in dt, we recover Eqs. l(T7|) 
and H lfcip . respectively. 

However, these more general forms for 9 it) and <j){t) 
lead to a more physical interpretation of our method. 
Since the dimensionless magnetization components along 
and perpendicular to B cS are M z = cos9, M x = 
sin 9 cos <b, and M„ = sin 6 sin <j>, we have 



M z it) 



M x (t) 
and 



taxih(t&nh- 1 [M z ito)]+a [ dr|B off (T)| 
V Jt 



= v/l-[M 2 (t)] 2 cos[0(t)], 



M y {t) = y/1- [M z (t)} 2 sm[0(t)}. 



(20) 
(21) 

(22) 



Independent of the coordinates, we must assure that 
(for the ith nanomagnet in the cth configuration) M 



changes its direction smoothly, in order to obtain a re- 
liable calculation for the overall M. Since each compo- 
nent of M cannot change dramatically, we must there- 
fore require 9 <C 2tt and <fi <§; 2ir. These restrictions 
then require us to set the time integration step width dt 
sufficiently small. If for example, j/a were on the or- 
der of 10+ 11 and \B cS \ were in the range 1CT 3 - 1CT 2 
T, we would require dt < 10~ n s. For sweep rate 
m 10" 2 T/s, where At = N t dt w 10" 4 s, N t must be 
on the order of 10 7 . Since we would need to recalculate 
the direction of the magnetization of each nanomagnet 
Nt times in each AB step, this would be a significant 
challenge with present day computers. 

One thing we can do to make our calculations feasi- 
ble for the sweep rates used in SMM studies is to set a 
extremely small, say a/7 < 10~ 10 , although such small 
a values have not been reported in experiments. Other- 
wise, to study much larger but perhaps more reasonable 
a values, we would have to use much faster sweep rates, 
as in KS.II1 
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